INTRODUCTION:

The metabolic disturbance of obesity is associated with premature mortality and it results from an imbalance between caloric intake and expenditure due to high-fat content in the diet^1,2. Data showed that approximately 108 million children and 608 million adults were considered obese in 2015. This number does not show a decreasing trend, as in 2018, in China alone, 85 million adults were considered obese³. Obesity is part of a cluster termed metabolic syndrome, which includes non-alcoholic fatty liver disease⁴ and chronic inflammation is a key feature of obesity⁵, and hypoxia has been identified as a potential risk factor for the chronic inflammation⁶. Obesity can be divided into monogenic, which results from a single genetic mutation, or polygenic obesity, which is a result of the interaction between the obesogenic environment and various genetic variants¹. The greatest challenge in any condition is not in the identification of the associated genes but in determining the molecular mechanisms by which the disease risk or the phenotypic expression can be reduced⁷. In our earlier work, it has been shown that the combination of different chemometric techniques of principal component analysis (PCA), random forest (RF), and genetic algorithm (GA) using multiarray data of high-fat diet fed mice, when combined with functional analysis, enabled the identification of a small subset of genes for diagnostic or therapeutic purposes⁸. The aim of this study is to extend the earlier work using data from liver samples instead of adipose tissue samples and to monitor if the approach identifies genes with new functional annotation capacity. The protein-protein interactions of the proteins identified by the three methods were to be compared, and important protein modules were ascertained. Finally, the aspect of regulatory control by microRNAs (miRNA) is to be extrapolated to human genes for identifying therapeutic or diagnostic targets.

MATERIAL AND METHODS:

Overview of Methods:

The methods’ workflow was similar to our earlier work⁸, but the data used was gene expression data from liver samples of high-fat diet fed mice instead of adipose samples. The genes were selected by principal component analysis, random forest, and genetic algorithm. The selected genes were evaluated for classification accuracy. The functional annotation of the selected genes and the protein-protein interactions were studied for each of the selected techniques.

Data

The dataset used in this study is the probe values for the liver samples from the earlier chosen dataset of GSE 39549. Microsoft Access was used to obtain the average expression of 15000 genes by mapping the probe sets of the genes in the raw data. The original data consisted of different time points for high-fat diet and normal diet treated animal livers, and these time points were pooled to compare a high-fat diet and a normal diet.

Software and packages

The selection of genes was by principal component analysis, random forest, and genetic algorithm. The free R package with prcomp library was used for PCA whereas RF was carried out using the library of randomForest. The selection of genes by the GA was done by using the software Matlab2019b. The libraries in R package of glm and e1071were used to study the ability of selected genes to classify the samples. The software Cytoscape 3.10.2 was used for network analysis and visualization and the analyses were carried out using an Intel CoreTM i9-12900 CPU@ 3.0 GHz with 64 GB RAM.

Selection of genes

The selection of fifty genes for the liver samples was similar to the earlier study with adipose samples⁸. The PCA was carried out using prcomp function in the R program, whereas the RF method using randomForest library had the mtry and ntree set to 120 and 1000, respectively. The GA was run using MATLAB with the same parameters as the earlier study⁽⁸⁾. The parameters used are the number of chromosomes of 100, ndims of 3, and the algorithm was iterated for 400 generations.

Machine learning for classification

The ability of the three sets of selected genes to differentiate the high-fat diet and control-diet samples were studied using six different machine learning techniques. However, only the first 30 ranked genes out of fifty were used for classification by k-nearest neighbors (kNN), and support vector machines.

Functional annotation clustering

The genes selected by PCA, RF, and GA were analyzed separately using the Functional Annotation tool in the Database for Annotation Visualization and Integrated Discovery or DAVID (https://david.ncifcrf.gov/). The biological meaning of the three sets of genes in terms of Gene Ontology (GO) functions was identified, especially those pertaining to biological processes, molecular functions, and cellular components related to lipid metabolism and obesity. The total 50 genes selected by each method were evaluated for functional annotation, and the similarity term overlap was set to 3. The similarity threshold was 0.5, and a p-value of less than 0.05 was used to obtain the optimal and statistically significant results. The enriched pathways of the genes, according to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, were also determined.

Protein-protein interactions network and hub genes

The protein-protein interactions in liver samples from high-fat diet treated animals were identified by submitting the three sets of genes to STRING: functional protein association networks (https://string-db.org/) database. The protein-protein interaction networks were determined using a confidence score of > 0.4, and the resulting networks were displayed by hiding the disconnected nodes. The chosen sources for the active interactions included “experiments, co-expression, neighborhood, gene fusion, co-expression, databases, and text-mining. The networks obtained were downloaded as tab-separated values format and further visualized in Cytoscape 3.72. A plug-in of Cytoscape, which is the Molecular Complex Detection or MCode, was used to obtain significant modules in each of the networks. The parameters of network scoring and clusters were node score cutoff =2, k-score = 2, max depth =100, cluster finding = haircut, and find clusters = in the whole network (9). The plug-in CytoHubba was used to identify hub genes as ranked by degree.

The associated miRNA-gene regulatory networks in human

The genes chosen by the three different multivariate methods were assessed for biological relevance in humans by studying the regulatory aspect of human miRNA on the human genes that were homologs of the mice genes identified by PCA, RF, and GA. The human protein-protein networks associated with the main mice protein modules were obtained using the STRINGIFY function of the STRING APP in Cytoscape. The CyTargetLinker was used to extend the human protein-protein network with the miRNA regulatory network from the database of miRTarbase and TargetScan v8⁽¹⁰⁾ to obtain the associated human miRNA-gene regulatory network. The regulatory aspects of genes selected by GA were done for the four-node linear human protein obtained by the STRINGIFY function.

RESULTS:

Genes differentially expressed in the liver from a high-fat diet and a normal diet

The PCA showed that principal component (PC) 1 contributed 27.3% of the overall variance. In contrast, PC 2 accounted for 17.6%, whereas other components played minor parts, with a total variance of 90% being reached by twenty-one PCs. The thirty genes chosen by PCA had the highest loadings or weightage for the first principal component. The first six of them were identified by their ENTREZ ID as Cyp2b9, Cidea, Gsta1, Thrsp, Gck, and Mgst3. The thirty most important genes for classification by RF were chosen based on the mean decrease in accuracy. The first six genes are Mgll, Sqle, Saa4, Tbc1d25, Prkcd, and Fh1. The out-of-bag (OOB) estimate of the error rate was 0% using number of trees of 500 and the number of variables tried at each split or mtry of 120. The choice of genes in GA was based on the ones with high loadings based on singular vector decomposition, and the first six important genes were Nog, Rptn, Pira11, Schip1, Palb2, and Rnf138rt1.

Gene ontology and pathway analyses

The set of genes identified by PCA identified 24 GO terms that had enrichment scores above one. Some of the terms related to metabolic processes or obesity were glutathione metabolic process, xenobiotic catabolic process, xenobiotic metabolic process, aromatase activity, iron-ion binding, steroid metabolic process, extracellular space, extracellular region, phosphatidylserine binding. The set of genes identified by RF identified 18 GO terms with enrichment scores greater than one. The terms related to metabolism were sterol biosynthetic process, lipid metabolic process, cholesterol metabolic process, cholesterol biosynthetic process, endoplasmic reticulum, and steroid metabolic process. The terms related to inflammatory processes were positive regulation of interleukin-1 beta production, positive regulation of inflammatory response, and the term of inflammatory response. The genes identified by GA had the least GO terms, with the number being 12, and none of the terms were directly associated with metabolic processes or obesity. Cell projection had 8 genes associated with it, whereas the term cell differentiation had 6 genes, whereas spermatogenesis had five terms. The KEGG pathway identified by genes chosen by PCA were chemical carcinogenesis-receptor activation, chemical carcinogenesis-DNA adducts, fluid shear stress and atherosclerosis, drug metabolism-cytochrome P450, glutathione metabolism, hepatocellular carcinoma, retinol metabolism as well as steroid hormone biosynthesis. The KEGG pathway identified by genes chosen by RF was steroid biosynthesis, whereas no KEGG pathway was identified by genes chosen by GA.

Protein-protein interaction networks, significant modules, and hub genes

The genes identified by PCA formed a PPI network consisting of 4 components, and the overall network consisted of 29 nodes and 44 edges. The network formed by genes selected by RF also consisted of 4 components, and it had 31 nodes and 41 edges (Fig 1). The MCODE plugin showed that the PPI network of PCA consisted of three modules, with the top one containing 6 nodes which were for the genes Cyp2b9, Gstm2, Gsta1, Cyp2b13, Mgst3, and Gsta4. The top module was related to the glutathione metabolic process. The second module had 4 nodes, which were Cidea, Cfd, Plin4, and Fitm1. The third module consisted of nodes Lcn2, Cxcl1 and Ccn2. The hub nodes in the protein network identified by PCA were Anxa5, Gsta1, and Cyp2b9, which had 7 degrees of connectivity. Two modules were identified for the RF PPI network, with the top one having 5 nodes which were Msmo1, Fdps, Sqle, Tm7sf2, and Stard4. This module was related to cholesterol metabolic as well as biosynthetic processes. The second module also had five nodes but had one edge less than the first module. The five nodes in the second module were Vcam1, Ciita, Casp1, Cd44, and Tlr2, and it was related to the inflammatory process. The top two hub nodes for RF PPI were Cd44 and Casp1, with 6 degrees of connectivity. The genes identified by GA formed a PPI network, which had 4 components. The overall number of nodes was 11, and the number of edges was 7. There were no protein modules identified by MCode, and no nodes had a degree of connectivity of three or more. The largest component in this network was a single linear chain of 4 nodes consisting of Palb2, Nme8, Eef1d, and Rpl36al.

Fig 1. Protein-protein interactions among genes chosen by random forest

Human regulatory networks

The top module of PCA, when converted to human PPI, showed that a similar human network consisted of the proteins GSTA1, GSTA4, GSTM2, and MGST3. The regulation of these proteins was by 27 human miRNAs. The gene MGST3 was regulated by 21 miRNAs and GSTA4 by 5 miRNAs. The gene GSTA1 was regulated by a human miRNA, hsa-miR-26b-5p. The gene GSTM2 was regulated by hsa-miR-181a-50. The top module of RF formed a network consisting of the 5 human genes linked with 165 human miRNAs. Most of the miRNAs were regulating MSMO1, with STARD4 having the second most miRNAs. The gene FDPS was regulated by three miRNAs, whereas SQLE was also regulated by three human miRNAs. The identification of human gene targets for regulation by human miRNAs for GA-selected genes was done for linear chain PPI of PALB2, NME8, EEF1D, and RPL36AL. The human gene EEF1D was regulated by 4 human miRNAs, whereas RPL36AL was regulated by three miRNAs. The gene NME8 was regulated by hsa-miR-484 and PALB2 was regulated by hsa-miR-221-3p.

DISCUSSION:

The use of bioinformatics to understand the biological meaning of subsets of gene variables identified with data mining techniques on gene expression data from adipose tissue has identified potential biomarkers for diagnosis and therapeutics of obesity⁸. This study compared the use of this technique on liver data and the results showed that a similar pattern of data reduction. The genes selected by RF were the best at classification of the different samples as well as finding biological meaning. The GA approach was the least capable of classification as well as identifying biological meaning. The GO terms for the PCA-identified genes were related to xenobiotics metabolism as well as steroid metabolism. The GO terms for the RF-identified genes were more related to inflammatory processes as well as lipid and cholesterol metabolic processes. These terms were not identified by genes selected by GA but it identified processes related to cell structures such as cilium and cell projection.

Another manner to identify key genes and proteins is the use of MCode to identify protein modules. This is because a module is a group of closely related proteins that act in concert to perform a specific biological function through PPI¹¹. The key modules of PCA and RF selected genes had key roles related to obesity. It has been reported that obesity-induced vascular dysfunction involves remodeling of visceral adipose tissue and increased inflammation. This process was accompanied by increased expression of Vcam-1¹², which is part of the second protein module of RF-selected variables. The other roles related to obesity can be attributed to genes such as Casp1 which is a hub node in the RF selected genes. Casp1 was found to be downregulated in human adipocytes when exposed to low oxygen or hypoxia¹³. It has been suggested that hub nodes have a high chance of taking part in essential biological processes⁽¹⁴⁾, and their deletion or dysregulation could be more deleterious to the organism¹⁵. The earlier work with adipose tissue samples showed that Vapa and Npc2 can act as crosstalk genes that link different processes⁸ but in this work the modules were apparently not connected. It was found that by combining the genes selected by PCA and RF the module related to glutathione metabolic process can be linked to cholesterol metabolic and biosynthetic process by the crosstalk genes of Abcg5, Abcg8, and Cyp3a11. The Abcg proteins play a major role in the maintenance of cellular cholesterol levels¹⁶. Overall, the genes identified by the three multivariate techniques in this study form part of the 241 upregulated and 91 downregulated genes in response to a high-fat diet¹⁷.

The identification of important genes for therapeutic or diagnostic purposes is confounded by the fact that miRNAs can regulate the expression and function of several genes¹⁸. The miRNAs have been reported to be differentially expressed in the adipose tissue of mice after a long-term high-fat diet¹⁹. This study showed that various miRNAs-gene interactions are present in humans for most of the genes identified by PCA, RF, or GA. Some genes that showed fewer interactions with human miRNAs are CIDEA, PLIN4, and NME8. The gene NME8 has been implicated in hepatocarcinoma²⁰, whereas Plin4 has been upregulated in the testis of high-fat treated rats²¹. MiRNAs are involved in regulating genes and they have been suggested as potential biomarkers as well as even therapeutic agents⁵. This study also shows the potential of extending the miRNA-regulatory model in animals to humans for predicting future research. Although GA found fewer genes that showed PPI, it is important to realize that this technique is suited for finding optimized solutions in a large search space. It could be suggested that the genes identified by it have novel mechanisms that are not yet discovered, or they interact less which could make them potential therapeutic targets. The advantage of using the approach of this study, where it highlights different sets of genes with biological function and regulatory mechanisms, should be tested on microarray or RNA-sequencing data from other diseases. The limitations of this study are that the number of genes studied is 15,000 and not the overall genes in the microarray data, the number of samples is still limited, and the regulatory aspects by miRNAs were not confirmed by experimental studies. Nevertheless, this study shows the potential of data mining to identify a set of genes that could be utilized for diagnostic or therapeutic purposes but with the caveat that further experimental work on human or animal samples are needed to validate the findings.

ABBREVIATIONS:

Abcg5: ATP-binding cassette sub-family G member 8*; Abcg8: ATP-binding cassette sub-family G member 5*; Anxa5: annexin A5*; Casp1: caspase 1*; Ccn2: cellular communication network factor 2*; Cd44: CD 44 antigen*; Cfd: complement factor D*; Cidea: cell death-inducing DNA fragmentation factor, alpha subunit-like effector A*; CIDEA: cell death-inducing DNA fragmentation factor, alpha subunit-like effector A**; Ciita: class II transactivator*; Cxcl1: C-X-C motif chemokine ligand 1*; Cyp2b13: cytochrome P450, family 2, subfamily b, polypeptide 13*; Cyp2b9: cytochrome P450, family 2, subfamily b, polypeptide 9*; Cyp3a11: cytochrome P450, family 3, subfamily a, polypeptide 11*; DAVID: Database for Annotation, Visualization and Integrated Discovery; EEF1D: eukaryotic translation elongation factor 1 delta##; Eef1d: eukaryotic translation elongation factor 1 delta*; EEF1D: eukaryotic translation elongation factor 1 delta**; Fdps: farnesyl diphosphate synthetase*; FDPS: farnesyl diphosphate synthetase*; Fh1: fumarate hydratase 1*; Fitm1: fat storage-inducing transmembrane protein 1*; GA: genetic algorithm; Gck: glucokinase*; GEO: gene expression omnibus; GO: Gene ontology; GSTA1: glutathione S-transferase, alpha 1##; Gsta1: glutathione S-transferase, alpha 1*; GSTA1: glutathione S-transferase, alpha 1**; GSTA4: glutathione S-transferase, alpha 4##; Gsta4: glutathione S-transferase, alpha 4*; GSTA4: glutathione S-transferase, alpha 4**; GSTM2: glutathione S-transferase, mu 2##; Gstm2: glutathione S-transferase, mu 2*; GSTM2: glutathione S-transferase, mu 2**; KEGG: Kyoto Encyclopaedia of Genes and Genomes; kNN: k-nearest neighbour; Lcn2: lipocalin 2*; Mgll: monoglyceride lipase gene*; MGST3: microsomal glutathione S-transferase 3##; Mgst3: microsomal glutathione S-transferase 3*; MGST3: microsomal glutathione S-transferase 3**; miRNA: microRNA; Msmo1: methylsterol monoxygenase 1*; MSMO1: methylsterol monoxygenase 1**; NME8: NME family member 8##; Nme8: NME family member 8*; NME8: NME family member 8**; Nog: noggin*; Npc2: Niemann-Pick type 2*; PALB2: partner and localizer of BRCA2##; Palb2: partner and localizer of BRCA2*; PALB2: partner and localizer of BRCA2**; PC: principal component; PCA: principal component analysis; Pira11: paired-Ig-like receptor A11*; Plin4: perilipin 4*; PLIN4: perilipin 4**; Prkcd: protein kinase C, delta*; RF: random forest; Rnf138rt1: ring finger protein 138, retrogene 1*; Rptn: repetin*; Saa4: serum amyloid A 4*; Schip1: schwannomin interacting protein 1*; Sqle: squalene epoxidase#; Sqle: Squalene epoxidase*; SQLE: Squalene epoxidase**; Stard4: StAR-related lipid transfer (START) domain containing 4*; STARD4: StAR-related lipid transfer (START) domain containing 4**; STRING: Search Tool for the Retrieval of Interacting Genes; SVM: support vector machines; Tbc1d25: TBC1 domain family, member 25*; Thrsp: thyroid hormone-responsive*; Tlr2: toll-like receptor 2*; Tm7sf2: transmembrane 7 superfamily member 2*; Vapa: vesicle-associated membrane protein, associated protein A*; Vcam1: vascular cell adhesion molecule 1* (*: mouse gene; **: human gene; #: mouse protein; ##: human protein)

CONFLICT OF INTEREST:

The authors declare no conflict of interest.

ACKNOWLEDGEMENT:

The data analyses approach of this project was carried out as part of project FRGS/1/2014/SKK01/ UNISZA/03/1. Dr Saravanan Dharmaraj acknowledges the financial backing of Ministry of Higher Education, Malaysia for the above research grant.

REFERENCES:

1. Loos RJF, Yeo GSH. The genetics of obesity: from discovery to biology. Nat Rev Genet. 2022; 23(2): 120–33.

2. Omer T. The causes of obesity: an in-depth review. Adv Obesity, Weight Manag Control. 2020;10(4): 90–4.

3. Wang L, Zhou B, Zhao Z, Yang L, Zhang M, Jiang Y, et al. Body-mass index and obesity in urban and rural China: findings from consecutive nationally representative surveys during 2004–18. Lancet [Internet]. 2021; 398(10294): 53–63. Available from: http://dx.doi.org/10.1016/S0140-6736(21)00798-4

4. Chen K, Shen Z, Gu W, Lyu Z, Qi X, Mu Y, et al. Prevalence of obesity and associated complications in China: A cross-sectional, real-world study in 15.8 million adults. Diabetes, Obes Metab. 2023; 25(11): 3390–9.

5. Ge Q, Brichard S, Yi X, Li Q. microRNAs as a new mechanism regulating adipose tissue inflammation in obesity and as a novel therapeutic strategy in the metabolic syndrome. J Immunol Res. 2014;2014:987285.

6. Ye J. Emerging role of adipose tissue hypoxia in obesity and insulin resistance. Int J Obes. 2009; 33(1): 54–66.

7. McCarthy MI, Abecasis GR, Cardon LR, Goldstein DB, Little J, Ioannidis JPA, et al. Genome-wide association studies for complex traits: consensus, uncertainty and challenges. Nat Rev Genet [Internet]. 2008 May; 9(5): 356–69. Available from: http://www.nature.com/articles/nrg2344

8. Dharmaraj S, Mahadeva Rao US, Simbak N. Advantages of functional analysis in comparison of different chemometric techniques for selecting obesity-related genes of adipose tissue from high-fat diet-fed mice. Biomed Res Ther. 2021; 8(3): 4277–85.

9. Zhang Y, Zheng Y, Fu Y, Wang C. Identification of biomarkers, pathways and potential therapeutic agents for white adipocyte insulin resistance using bioinformatics analysis. Adipocyte [Internet]. 2019; 8(1): 318–29. Available from: https://doi.org/10.1080/21623945.2019.1649578

10. Chai Y, Tan F, Ye S, Liu F, Fan Q. Identification of core genes and prediction of miRNAs associated with osteoporosis using a bioinformatics approach. Oncol Lett [Internet]. 2018 Sep 26; 17(1): 468–81. Available from: http://www.spandidos-publications.com/10.3892/ol.2018.9508

11. Lin CY, Lee TL, Chiu YY, Lin YW, Lo YS, Lin CT, et al. Module organization and variance in protein-protein interaction networks. Sci Rep. 2015;5:1–12.

12. Yao L, Bhatta A, Xu Z, Chen J, Toque HA, Chen Y, et al. Obesity-induced vascular inflammation involves elevated arginase activity. Am J Physiol - Regul Integr Comp Physiol. 2017; 313(5): R560–71.

13. Mazzatti D, Lim F-L, O’Hara A, Wood IS, Trayhurn P. A microarray analysis of the hypoxia-induced modulation of gene expression in human adipocytes. Arch Physiol Biochem [Internet]. 2012; 118(3): 112–20. Available from: http://www.tandfonline.com/doi/full/10.3109/13813455.2012.654611

14. Ouyang S, Tang R, Liu Z, Ma F, Li Y, Wu J. Characterization and predicted role of microRNA expression profiles associated with early childhood obesity. Mol Med Rep. 2017; 16(4): 3799–806.

15. He X, Zhang J. Why do hubs tend to be essential in protein networks? PLoS Genet. 2006; 2(6): 0826–34.

16. Sun Y, Wang J, Long T, Qi X, Donnelly L, Elghobashi-Meinhardt N, et al. Molecular basis of cholesterol efflux via ABCG subfamily transporters. Proc Natl Acad Sci U S A. 2021; 118(34): 1–10.

17. Do GM, Oh HY, Kwon EY, Cho YY, Shin SK, Park HJ, et al. Long-term adaptation of global transcription and metabolism in the liver of high-fat diet-fed C57BL/6J mice. Mol Nutr Food Res. 2011; 55(SUPPL. 2): 173–85.

18. Kakati T, Bhattacharyya DK, Barah P, Kalita JK. Comparison of Methods for Di ff erential Co-expression Analysis for Disease Biomarker Prediction. Comput Biol Med [Internet]. 2019; 113(May): 103380. Available from: https://doi.org/10.1016/j.compbiomed.2019.103380

19. Chartoumpekis D V., Zaravinos A, Ziros PG, Iskrenova RP, Psyrogiannis AI, Kyriazopoulou VE, et al. Differential expression of microRNAs in adipose tissue after long-term high-fat diet-induced obesity in mice. PLoS One. 2012; 7(4).

20. Liu G-M, Zeng H-D, Zhang C-Y, Xu J-W. Key genes associated with diabetes mellitus and hepatocellular carcinoma. Pathol Res Pract [Internet]. 2019; 215(11): 152510. Available from: https://doi.org/10.1016/j.prp.2019.152510

21. Yang X-Y, Gu Y-J, An T, Pan Y-Y, Mo F-F, Miao J-N, et al. Proteomics Analysis of Testis of Rats Fed a High-Fat Diet. Cell Physiol Biochem [Internet]. 2018; 47: 378–89. Available from: https://www.karger.com/Article/FullText/486732

Received on 27.06.2024 Modified on 30.07.2024

Research J. Pharm. and Tech 2024; 17(8):4043-4048.

DOI: 10.52711/0974-360X.2024.00627